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This  article  concerns  the  effects  of  non-Fickian  water  diffusion  in  £ber- 
reinforced  polymeric  composites.  The  departure  from  classical  diffusion 
is  attributed  to  the  time-dependent  response  of  the  polymer ,  akin  to  vis¬ 
coelastic  mechanical  response. 

A  formulation  is  provided  to  evaluate  the  coefficient  of  moisture  diffu¬ 
sion  and  moisture  profiles  within  the  composite  for  the  non-Fickian  case. 
In  addition,  it  is  demonstrated  that  computed  magnitudes  of  residual 
hygro-thermal  stresses  may  differ  by  about  25%  from  predictions  based 
upon  classical  diffusion. 


1.  Introduction 

The  diffusion  of  water  in  polymers  and  fiber-reinforced  polymeric 
composites  has  been  studied  extensively  for  almost  a  century  and  a 
vast  body  of  literature  deals  with  that  subject.  A  survey,  primarily 
concerned  with  fiber-reinforced  polymeric  composites,  was  published 
recently  (Weitsman,  1991). 
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The  most  prominent  and  common  formulation'  of  the  diffusion 
process  employs  the  well  known  Fick’s  law  (Fick,  1855),  whose  one 
dimensional  version  reads 

dC  d2C 

‘>0,~L<z<L.  (1) 

The  field  equation  (1)  is  accompanied  by  initial  and  boundary 
conditions 

C(z,  0)  =  Ci(z),  —L  <  z  <  L,  (2) 

and,  say 

C(±L,t)  =  CL(t),  t>  0.  (3) 

In  equations  (1),  (2)  and  (3),  z  and  t  are  spatial  coordinate  and 
time,  respectively,  C  =  C(z,  t )  is  moisture  content  and  D  denotes  the 
coefficient  of  moisture  diffusion.  In  many  cases  the  coefficient  D  is 
assumed  constant,  whereby  the  diffusion  process  is  Unear.  The  above 
assumption  will  be  employed  in  the  present  work. 

Consider,  in  addition,  the  special  circumstance  of  an  initially  dry 
plate  subjected  to  constant  boundary  conditions,  namely, 

C(z,0)  =  0,  —L  <  z  <  L,  (2a) 


C(±L,t)  =  t>  0.  (3a) 

where,  H(t )  is  the  Heaviside  step  function. 

Denote  the  solution  to  the  unit  step  input  ( i.e .,  Co  =  1  in  equation 
(3a))  by  Cjj(z ,  t).  The  well  known  result,  given  in  Crank  (1975),  reads: 


Ch(zj) -  j>ir  b 


+  erfc 


2 n  -f- 1  -f -  z/ L 


Integration  across  the  thickness  provides  the  weight  gain  Mh  (f),  which 
corresponds  to  Cn(z,t), 


Mn(t)  =  4 Ly/t*  • 


(5) 


2 


where 


r." 

L2  ’ 


(6) 


erfc(^)  is  the  complementary  error  function  and  ierfc(x)  is  its  integral: 

ierfc(ar)  =  /  erfc  (£)cf£.  (7) 

Jo 


In  view  of  linearity,  the  solution  for  Co  /=■  1,  as  well  as  for 
C(±L ,  t)  =  g(t),  can  be  expressed  by  means  of  Ch(z,  t)  and  Mn(t)  in 
a  straightforward  manner  (Crank,  1975). 

Since  data  on  moisture  distribution  are  scarce  and  difficult  to  gen¬ 
erate,  the  most  readily  available  experimental  information  accounts 
for  total  weight  gain.  Before  comparing  weight  gain  data  with  model 
predictions,  note  the  characteristic  features  of  when  plotted 

vs.  as  shown  in  Figure  1.  Accordingly,  the  value  of  the  initial 
slope  is  4 L/yft,  departures  from  linearity  by  ±1%  do  not  occur  until 
y/t*  «  0.557,  Mn(t*)  «  0.62,  and  saturation,  to  within  ±1%,  occurs 
at  y/t*  =  1.32.  It  should  be  recognized  that  typical  weight  gain  data 
exhibit  scatter  of  at  least  ±1%. 

In  addition,  it  is  possible  to  infer  the  value  of  the  diffusion  coeffi¬ 
cient  D  from  (Shen  and  Springer,  1981),  we  have 


D  = 


~lim 
16  t-o 


t 


(8) 


In  this  manner,  knowledge  of  enables  the  evaluation  of 

C(z,t),  through  retracing  expression  (4)  from  equation  (5),  and  the 
computation  of  residual  hygral  stresses  by  means  of  well  established 
analyses  (Tsai  and  Hahn,  1980,  Harper  and  Weitsman  1985,  Weits- 
man,  1991). 

Several  uncertainties  arise  when  moisture  weight  gain  data  do 
not  correspond  to  the  strictures  exhibited  in  Figure  1,  implying  that 
the  premises  which  led  to  Mn(t)  are  not  met  by  the  polymeric  com¬ 
posite  material  at  hand.  Such  circumstances  occur  very  frequently, 
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with  a  typical  example  exhibited  in  Figure  2  for  the  case  of  Fiberite 
T300/1014  composite  (Blikstad,  et  al ,  1988).  To  accentuate  the  de¬ 
partures  from  “classical”  predictions,  which  corresponds  to  equation 
(5),  the  data  in  Figure  2  are  bracketed  by  two  curves  which  conform 
to  the  format  of  Figure  1  (multiplied  by  two  distinct  constants  Cq  to 
provide  best  fits  with  data  at  short  and  long  times,  respectively).  For 
future  reference,  these  curves  are  denoted  by  “Upper  Fickian”  and 
“Lower  Fickian”,  respectively. 

There  are  several  plausible  explanations  for  the  causes  for  depar¬ 
ture  from  a  linear  diffusion  process  (Weitsman,  1991).  Of  these,  one 
rationale  is  motivated  by  considerations  akin  to  the  well  known  time 
dependent,  viscoelastic  response  of  polymers.  Accordingly,  the  very 
same  Gibbs  free  energy  which  gives  rise  to  time-dependent  mechani¬ 
cal  response,  predicts  a  diffusion  process  with  time-dependent  bound¬ 
ary  conditions  even  under  exposure  to  constant  ambient  environment 
(Weitsman,  1990). 

This  consideration  will  be  employed  in  the  present  work,  namely, 
we  shall  retain  the  field  equation  (1)  and  initial  condition(2a)  but  will 
consider  the  boundary  condition 

C(±L,  *)  =  /(*)  (9) 

instead  of  (3a),  even  though  the  ambient  condition  remains  constant 

The  following  issues  will  be  addressed  in  the  present  work: 

(a)  Express  the  moisture  distribution,  C(z,t ),  when  moisture  up¬ 
take  M(t)  does  not  conform  to  expression  (5). 

(b)  Establish  a  method  to  determine  the  value  of  the  coefficient 
of  moisture  diffusion,  D,  when  expression  (5)  —  and  thereby  also 
equation  (8)  —  no  longer  hold. 

(c)  Evaluate  the  effects  of  foregoing  distribution  C(z ,  t )  on  residual 
stresses  in  composite  laminates. 
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2.  Diffusion  with  Boundary  Conditions 
of  Viscoelastic  Type 


It  has  been  suggested  ( e.g .,  Long  and  Richman,  1960)  that  depar¬ 
ture  from  “classical”  diffusion  may  be  explained  by  replacing  bound¬ 
ary  condition  (3a)  with 

C(±L,  t )  =  [C0  +  Ci  (l  -  e-*)]  H(t).  (10a) 

More  recently,  it  has  been  shown  (Weitsman,  1990)  that  for  vis¬ 
coelastic  materials,  both  creep  compliance  S(t )  and  chemical  potential 
p(t)  are  expressible  in  Prony  series  forms,  namely, 


/*(*)  =  Mo  +  J2  Mn  (l  -  e  pnt) 

n— 1 


This  observation  suggests  that,  in  view  of  the  time-dependent  re¬ 
sponse  of  the  polymer,  one  should  consider  the  boundary  condition 


C(±L,t)  = 


Co  +  EC7„(l-e-^) 

n=l 


H{t) 


(10b) 


of  which  (10a)  accounts  for  the  first  term  in  the  series.  The  above 
expression  implies  that  equilibrium  between  the  moisture  content  just 
inside  the  material  and  the  chemical  potential  of  the  external  vapor 
is  established  over  an  extended  time  —  rather  than  instantaneously. 

The  solution  to  equation  (1),  with  initial  condition  (2a)  and 
boundary  condition 


C(±L,t)=  (l  -  e"*1)  H(t) 

is  well  known  (Crank,  1975).  For  future  reference,  it  will  be  denoted 
by  C(z,t;0).  We  have 


C(z,t-,p )  =  1  -exp(-pt) 


COS  Zy/$5 

COsLy/fiD 


1 6/3L2  ^2,  (— l)nexp{— [(2n  +  l)7r/2]2t*}  (2n  +  l)nz 


“  (2 n  +  1)[4 PL2  -  Dir(2n  +  l)2] 


cos 


2  L 


,(11) 


5 


Upon  integration  across  the  thickness,  the  total  weight  gain  corre¬ 
sponding  to  C(z ,  t ;  ft)  is 

M{t\p)  =  2L-  |l  —  exp(-/9<)^^jtan 

8  y,  exp{— [(2n  +  l)7r/2]2t*} 

"^n=o  (2n  +  1)2{1  -  (2n  +  l)a[0,rV(4/?Za)]} 

Consequently,  the  distribution  and  total  weight  gain  due  to 
boundary  condition  (10b)  are 

C(M)  =  C0CH(M)  +  £  CnC(z,t\0n),  (13) 

n=l 

and 

M(t)  =  c„MH(f)  +  E  (14) 

n=l 

The  latter  expression  should  be  correlated  with  experimental 
data. 

3.  Data  Fitting 

Since  data  contain  experimental  error  and  statistical  scatter,  it  is 
advisable  to  smoothen  the  weight  gain  data,  Mexp(t),  before  affecting 
a  match  with  expression  (14). 

Noting  that  for  early  times  M#(i)  is  proportional  to  y/t,  we  choose 
to  fit  the  experimental  data  as  follows, 

Mfit(<)  =  EV/2.  (15) 

j=i 

The  coefficients  Aj  are  determined  by  minimizing  the  square  error, 
namely, 

J  I  2 

gj-jf'"'"  JUO-EV'1  <»  =  0,  (t  =  1,2,---, J)  (16) 
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In  equation  (16),  tmax  denotes  the  duration  of  the  moisture  uptake 
experiment. 

With  Aj  determined  through  expression  (16)  it  is  now  necessary 
to  express  Mfit  in  the  form  (14).  This  task  is  accomplished  in  two 
steps. 

Consider  first  the  step-wise  increasing  boundary  condition 
C(±L,t)  =  C0H(t)  +  £  C?H(t  -  tk), 

k= 1 

to  which  corresponds  the  total  weight  gain 

M{t)  =  C0M„(t)  +  £)  C”M„(t  -  tk)H(t  -  tk).  (17) 

k=l 

The  quantities  Co,  C £  and  tk  (k  =  1,2 ,•••,#)  in  equation  (17) 
can  be  selected  in  a  manner  which  yields  M(t)  —  <  SCT  at  all 

times  f,  with  <5cr  a  prescribed  tolerance.  It  is  advisable  to  select  SCT  to 
be  smaller  than  a  typical  discrepancy  between  Mexp(<)  and 

The  proposed  procedure  is  to  select  Co  in  equation  (17)  which,  in 
view  of  equation  (5),  would  yield  M(t)  =  Mfit  at  the  first  time  step, 
i.e., 

_  Mfit(A<)  j  7T 
°“  4  V  DAt' 

where  At  is  the  selected  time  step. 

Then,  retain  M(t)  =  CoMjj(t)  imtil  such  time  t  =  t\  when 
Affit(fi)  —  CoMf{(ti)  =  6  «  SCT.  In  view  of  equation  (5),  this  dis¬ 
crepancy  at  t  =  t\  is  overcome  by  introducing  a  step  increment  in  the 
boundary  condition  at  some  earlier  time  t\  of  magnitude 


The  combined  effect  of  CoM//(f)-f-Cf  Mn(t  —  ti)  is  now  compared 
with  Mexp(t)  until  such  time  t2  when  they  differ  by  <5cr,  at  which  stage 
another  incremental  step  is  added  to  the  boundary  condition  at  an 
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earlier  time  <2  of  magnitude  C$  analogous  to  Cf .  The  procedure  is 
repeated  until  the  entire  range  of  Mexp(<)  is  covered.  The  amplitudes 
of  the  incremental  steps  are  given  by 


S  I  7?  ~~ 

4  \j~D(ik  -tVy 


(18) 


In  our  computations,  we  found  it  is  expedient  to  select  4  to  co¬ 
incide  with  ?jfc_i  (with  t\  =  0). 

With  known  C\ f  and  tk  (k  =  1, 2,  •  •  ♦ ,  A),  it  is  possible  to  convert 
the  step-wise  incrementing  boundary  condition 


fc=i 


(19) 


to  the  continuous  Prony  series  form 


c(±i,t)  =  £  c„  (1  -  e’e"') 

n=l 


(20) 


Employing  a  least  square  fit,  we  have 


(21) 


Equation  (21)  results  in  an  iV  x  iV  system  of  linear  algebraic 
equations  in  Cn 


a.jCj  = 


(22) 


where 


aij  —  ^max 


Pi 


Pi 


Pi  +  Pj 


K 

*.-E 

*=i 


cF 


(23) 

(24) 


It  should  be  noted  that  the  above  procedure  provides  some  lati¬ 
tude  in  the  selection  of  4  (k  =  1, 2,  •  •  • ,  A)  and  /3j  ( j  =  1,2,---,  N ). 
Although  no  hard  and  fast  rule  seems  to  be  available,  it  is  advisable 
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to  select  (3j  which  cover  the  spectrum  of  experimental  time,  preferably 
two  —  but  at  least  one  —  values  of  /3  per  time  decade. 

It  may  appear  that  the  employment  of  the  intermediate  series 
(19)  is  redundant  and  that  the  values  of  Cn  in  equation  (20)  can  be 
obtained  directly  through  a  least  square  fit  of 


d 

dC 


2 

r  /  Mfit(t)  -  C0M„(t)  -  £  CnM(t;  /?„)  , 

i  J0  n=l 


dtt  =  0, 


(25) 


which  yields  an  N  x  N  system  of  linear  algebraic  equations 


a  &  =  k 


(26) 


to  determine  Cj. 

In  this  circumstance  we  have 


and 


ft  mu  A  A 

a0  =  /  M{t-,  pj)dt, 

Jo 

i,  =  [jw*(<)  -  c„mh(<)]  mu  0,)dt. 

JO 


(27) 


(28) 


Unfortunately,  the  numerical  evaluation  of  and  6,  involves  com¬ 
pounded  effects  of  truncation  errors  and,  most  critically,  yields  an  ill- 
conditioned  matrix  d,j.  The  latter  difficulty  arises  from  the  fact  that 
rows  and  columns  in  the  a,j  matrix  consist  of  elements  of  very  close 
magnitudes. 

In  many  circumstances  moisture  uptake  data  do  not  seem  to  ap¬ 
proach  an  equilibrium  value,  regardless  of  the  duration  of  exposure, 
and  M(t)  tends  to  increase  according  to  M{t)  ~  Ktp  as  t  »  1.  It 
can  be  shown  that  this  circumstance  is  commensurate  with  a  time 
dependent  boundary  condition 


C(±L,t)  =  C0rH(t). 


(29) 


To  prove  this,  let  C(z ,  s )  denote  the  Laplace  transform  of  C(z,t), 
then,  since 

r(p+i) 


C(±L,s)  =  Co- 


sp+l 
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it  can  be  shown  that 


C(z,s)  = 


C0r(p+1)  coshgz 
sp+1  cosh  qL 


(30) 


where  q  =  \Js/D. 

Consequently,  the  Laplace  transform,  M(s),  of  the  total  weight 
gain  M(t)  is 

2C0r(p  + 1)  tanhgL 


M{s)  = 


(31) 


sp+1  q 

Although  it  seems  impossible  to  express  M(t)  analytically1,  its 
asymptotic  value  for  t  »  1  is  readily  obtainable  from 


which  yields 


limM(s)  =  2L-°T^+}.\ 
s— *0  V  '  sp+1 


lim  M(t)  =  2 LC0tp. 

t-fOO 


(32) 


4.  Evaluation  of  the  Coefficient  of  Moisture  Diffusion 

In  all  previous  computations,  it  was  implicitly  assumed  that  the 
value  of  the  coefficient  of  moisture  diffusion  D  is  known.  However, 


1  An  exception  occurs  for  the  important  circumstance  p  =  1/2,  which  corresponds 
to  M(t)  ~  Ay/t  as  t  —*  oo.  In  this  case  both  C(z,t)  and  M(t)  have  the  following 
analytical  expressions  (see  Appendix) 


C(z,t )  =  CQy/t 


(~l)n+1  „  ( (2n  +  1) 


i  — - —  y*  -i— ^ — f 

7T2V^^0(2n  +  l)2 


M(t)  =  2LC0Vt 


L  2  v< )\- 


+  l)nz 


where  F(()  denotes  Dawson’s  integral  (Olver,  1974)  defined  by 


F(0  =  e-«2  /  e“adu. 
J  0 
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when  experimental  weight  gain  data  do  not  fit  the  format  shown  in 
Figure  1,  expression  (8)  does  not  apply  and  D  is  unknown. 

Assume  as  before  that  departure  from  classical  diffusion  are  at¬ 
tributable  strictly  to  the  time-dependence  of  the  boundary  condition, 
namely, 

C(±L,t)  =  C0[l  +  f(t)]H(t)  (33) 

with  /( 0)  =  0. 

Therefore, 


M(t)  =  Co  +  jf  MH(t  -  . 

Let  g(s )  denote  the  Carson  transform  of  g{t ),  namely, 


then 


whence 


yoo 

g(s)  =  s  e~atg{t)dt , 

Jo 


M(s)  =  Co  [M*(s)  +  M„(s)f(s)\  , 


^  CoMh(s) 


Note  that  at  this  stage,  Co,  f(t)  (whereby  also  f(s ))  as  well  as  D 
are  unknown.  However,  Mh{s )  is  known  analytically  (Crank,  1975), 
namely, 

M„(s)  =  2^tanh]fjjL.  (36) 

In  addition,  M(s),  which  is  the  Carson  transform  of  the  weight 
gain  data,  can  be  computed  numerically. 

By  hypothesis,  f(t )  is  independent  of  the  sample  thickness  L. 
Therefore,  consider  two  sets  of  weight  gain  data  M (t;  L\)  and  M(t\ L2) 
associated  with  thicknesses  L\  and  L2,  respectively.  Equation  (35) 
yields 

=  Mh^U) 

M(s;L2;D)  Mh(s;L2)'  K  ’ 
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Note  that  equation  (37)  contains  the  unknown  D  alone  and  does 
not  include  f(s)  and  Co- 

Consider  a  data  fit  according  to  equation  (15),  whose  Carson 
transform  is  given  by 

*  W  =  E  4r4  + 1).  (38) 

j- 1  s  z 

Consequently,  it  is  possible  to  compute  the  ratio 

pk  =  p(Sk )  = 

P  M(sk;L*,D) 

at  distinct  values  of  transform  parameter  s  =  Sk  (k  =  1,2,  •  •  ■ ,  K)  and 
express  the  ratio 


rk(D)  =  r(sk,D)  = 


tanh  \jskL\i 
tanh  \JskL\i 


which  corresponds  to  the  right  hand  side  of  equation  (37). 

The  value  of  D  can  be  determined  from  the  best  least  square  fit 
obtained  from 

d  JL 


namely, 


where, 


3nEk-rt(D)f  =  0, 


£,[Pk  -  rk(D))r'k(D)  =  0, 


r'k(D )  = 

-Licsch2(LiyJsk/ D)tsnh(Li\/sk/ D)]. 


Denote 


hra)* 


9k{D)  = 


(41) 


=  L2cth(L1^Sk/D)sech2(L2\/sk/D)  - 
-Licsch2(Liyfsk/D)  tanh(Zi  y/sk/D), 

then,  since  D  ^  0,  eqn.(40)  is  equivalent  to 

I >*  ~  rk(D)]gk(D)  =  0.  (42) 

jt=i 

This  equation  can  be  solved  numerically  for  D.  In  the  sequel, 
Newton’s  Iteration  Method  will  be  employed. 

Newton’s  Iteration  Method  starts  with  a  suitably  chosen  initial 
value  P(°)  and  employs  iteration  until  attaining  convergence  within  a 
prescribed  tolerance.  In  the  present  case,  the  nth  iterative  value,  D^n\ 
is  related  to  the  (n  —  l)st  value  as  follows: 


K 


Ek-#-1')]#-1’) 

£)(n)  _  £)(n-l)  _  *=1 


'  lD=D<n~1) 

The  denominator  in  equation  (43)  can  be  expressed  as 


4 


(43) 


{K  s  t 

X>*  -  =  E  (b*  -  r*(DM(D)  +  H^(D)} 


where 


gk(D )  =  j  •  {cth(Liy/sk/D)tanh(L2yJsk/D) 


L\ 


[sinh2(Ii^/^) 

cosh2 (L2y/sk/D)\  sinh2(LiyJsk/D)  cosh2(L2y/sk/D) } 


L2 


L\L2 


5.  Moisture  Effects  on  Residual  Stresses 

To  demonstrate  the  significance  of  non-Fickian  diffusion  on  resid¬ 
ual  stresses  in  composite  materials,  consider  the  case  of  a  [O^/QO^js 
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symmetric  lay-up.  The  basic  stress-strain  relation  for  the  cross-ply 
laminate  are  (Harper  and  Weitsman,  1985): 

For  the  0°  layers, 

=  QL[t°x(t)  + znx(t)  -  aLAT- pLCe(z,t)] 

+  zKy(t)  -  oct AT  -  pTCe(z,  <)],  (44) 

cry(z,t)  =  QLT[t°x(t)  +  ZKx(t)  -  aLAT  -  PLCe(z,t)] 

+Qr[f°(<)  +  zKy(t)  -  qtAT  -  PjC^z^t )],  (45) 

where  subscripts  L  and  T  denote  the  the  directions  along  and  per¬ 
pendicular  to  the  fibers  respectively,  x  and  y  denote  the  length  and 
width  directions  of  the  laminate  respectively,  ck’s  are  the  thermal  ex¬ 
pansion  coefficients  and  /?’ s  the  moisture  swollen  coefficient,  AT  is  the 
temperature  variation,  Ce(z,t)  the  effective  moisture  concentration. 

For  the  90°  layers,  the  stresses  can  be  obtained  by  interchanging 
subscripts  L  and  T. 

In  the  above  equations,  e°  and  k  are  determined  in  terms  of  the 
external  applied  loads.  In  the  absence  of  such  loads  we  have 


Nx{t)  = 

/  <7x(z,t)dz  =  0, 

(46) 

# 

II 

rL 

/  ay(z,t)dz  =  0, 

'  ~  L 

(47) 

II 

5 

fL 

1  ax(z,t)zdz  =  0, 

“X/ 

(48) 

ll 

5 

fL 

1  ay(z,t)zdz  =  0. 

~L 

(49) 

Substitution  of  the  appropriate  stress  expressions  into  eqns.(46)-(49) 
yields,  after  some  manipulations,  the  following  results: 

4(0  =  £{( QL  +  Qr)[hPAT(t)  +  RG(t)  +  SF(t )] 

-QLrL[hPAT(t)  +  RF(t )  +  SG(t)]} ,  (50) 

*5(0  =  |  {Ql  +  QT)[hPAT(t)  +  RF(t)  +  SG(t)] 

-  QLTL{hPAT(t )  +  RG(t)  +  SF(t)}}  ,  (51) 
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where 


«*(<)  =  Ky(t)  =  0, 


(52) 


A  =  h\QL  +  Qt)(Ql  +  Qt )  -  ( QltL)\  F{t)  =  /*  Ce(*,  <)dz, 

./o 

/L 

G(t )  =  /  Ce(2,  t)dz,  P  =  +  Qltolt  +  Qlt&l  +  Qt<Xt, 

J  h 

R  —  QlPl  +  QltPt,  S  —  QltPL  +  QtPT,  Ce{z,t)  =  C{z,t)  -  Cth, 

Cth  denotes  the  threshold  moisture  concentration  below  which  no 
swelling  effect  is  observed  and  h  is  the  common  thicknesses  of  90° 
and  0°  laminae  group,  while  —  as  before  —  2 L  is  the  total  thickness 
of  the  laminate. 


6.  Computational  Results 

The  moisture  weight  gain  data  of  Figure  2  were  fitted  by  foregoing 
numerical  scheme  with  a  smooth  curve  Mfa(t)  according  to  equation 
(15)  and  subsequently  with  M(t)  according  to  equation  (14)  and  (5). 

Those  computations  yielded  a  six-term  fit  with  A\  =  0.23281, 
A2  =  0.010999,  Az  =  -7.2414  x  10“3,  A4  =  1.5040  x  10~4,  A5  = 
—1.3487  x  10"6,  Aq  =  4.4508  x  10~9  for  Mfit  and  the  resulting  curve 
is  shown  by  the  solid  line  in  Figure  4.  The  smooth  curve  of  Mfit  was 
subsequently  represented  by  M (t)  of  equation  (14)  with  6  =  O.OlMoo- 
This  representation  is  shown  in  Figure  4. 

Afterward,  the  data  were  represented  by  M(t)  as  expressed  by 
equation  (5)  with  Co  =  0.7  and  N  =  7,  C\  =  0.098605,  C2  =  0.021929, 
C3  =  -0.083097,  C4  =  0.27450,  C5  =  -0.55194,  C6  =  1.7834,  C7  = 
-1.5819  and  ft  =  1/10,  ft  =  1/50,  ft  =  1/100,  ft  =  1/500,  ft  = 
1/1000,  ft  =  1/5000,  ft  =  1/10000  (/3  in  hours"1). 

In  all  above  computation,  we  took  D  =  3.2  x  10-4  mm2/hour. 

The  resulting  curves  for  M(t)  are  indistinguishable  from  Mat(t) 
and  coincide  with  the  solid  line  shown  in  Figure  2. 
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For  purpose  of  comparison,  the  data  of  Figure  2  were  bracketed 
by  two  curves  of  the  form  CoMjj(t).  The  value  of  Co  =  1.12  x  10“2, 
with  D  =  3.2  x  10-4  mm2 /hour,  provide  a  “Lower  Fickian  Fit”,  while 
Co  =  1.32  x  10“2,  with  D  =  1.7  x  10-4  mm2/hour,  gave  an  “Upper 
Fickian  Fit”  to  the  data.  These  curves  are  shown  by  dashed  lines  in 
Figure  2. 

Employing  the  aforementioned  values  of  C,-  and  fa  (i  =  1, 2,  •  •  • ,  8) 
and  Co  =  0.7  x  10~2,  the  distributions  C(z,t)  were  computed  at  times 
t  —  500, 3600  and  8000  hours.  Results  are  shown  in  Figures  5  through 
7,  where  moisture  profiles  associated  with  both  “Upper”  and  “Lower” 
Fickian  data  fits  are  shown  for  comparison. 

Residual  hygro-thermal  stresses  were  evaluated  according  to  equa¬ 
tions  (44)  and  (45).  The  computations  employed  the  manufacturer’s 
data  (ICI  data  sheet)  which  gave  h  —  0.127mm,  El  =  148  GPa 
and  Et  =  11  GPa.  In  the  absence  of  further  details,  we  assumed  a 
Poisson  ratio  vLt  =  0.28.  The  above  quantities  were  converted  to 
the  stiffness  utilized  in  equations  (44)  and  (45)  as  follows  (Tsai  and 


Hahn,  1980):  vTL  =  vLT-~,  QL  —  - - - - 

El  1  —  vlt^tl 


>  Qt  = 


Et 


1  —  vlt^tl 


Qlt  = 


ultEt 

1  —  VLTVTL 


Furthermore,  we  assumed  hygrothermal  properties  similar  to 
those  of  AS4/3502  (Harper  and  Weitsman,  1985).  We  thus  took 
Pl  =  0,  Pt  =  3.24  x  10"3/1%  moisture  weight  gain,  Cth  =  0.1%, 
AT  =  —  180°C,  aL  =  -0.02  x  10-6/°C,  aT  =  27.5  x  10-6/°C. 

Results  of  the  hygro-thermal  residual  stresses  are  shown  in  Figures 


8-13. 


It  can  be  noted  from  the  above  results  that,  at  early  stages,  al¬ 
though  all  the  predicted  moisture  distributions  correspond  to  the  same 
total  moisture  weight  gain,  they  yield  distinct  stress  distributions.  At 
later  stages,  as  the  moisture  distribution  tends  be  uniform,  the  pre- 
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dieted  residual  stresses  are  basically  proportional  to  the  predicted  to¬ 
tal  moisture  weight  gain.  At  some  locations,  the  Fickian  predictions 
may  differ  by  as  much  as  20-30%  from  the  non-Fickian  values. 

Furthermore,  at  the  early  stages,  the  stresses  are  closer  to  the 
lower  Fickian  fit  since  the  boundary  moisture  saturation  level  is  closer 
to  the  lower  bound,  while  at  later  stages,  as  the  moisture  saturation 
level  increases,  the  stresses  are  closer  to  those  which  correspond  to 
the  upper  bound.  In  any  case,  it  is  worth  noting  that  the  lower  and 
upper  Fickian  fits  on  the  moisture  distribution  give  lower  and  upper 
bounds  of  the  resulting  hygro- thermal  stress  predictions. 

7.  Conclusions 

In  principle,  the  diffusion  process  of  water  in  polymers  and  poly¬ 
meric  composites  may  depart  from  the  idealizations  inherent  in  the 
classical  formulation  of  Fick.  Although  there  are  abundant  reasons  for 
that  departure,  the  current  work  focused  on  correlating  non-Fickian 
weight  gain  data  with  a  time-dependent  boundary  condition,  as  mo¬ 
tivated  by  the  viscoelastic  response  of  polymers. 

It  is  shown  that  non-Fickian  aspects  of  the  moisture  absorption 
process  have  significant  effects  on  the  resulting  hygro- thermal  stresses 
in  composite  materials.  The  current  work  reveals  the  sensitivity  of  the 
residual  hygro-thermal  stresses  to  a  boundary  condition  of  a  viscoelas¬ 
tic  type,  exhibiting  discrepancies  of  about  20-30%. 

For  a  diffusion  process  whose  lower  and  upper  Fickian  fits  can  be 
found,  is  it  safe  to  say  that  the  residual  stresses  which  correspond  to 
those  two  cases  provide  bounds  for  the  non-Fickian  case,  resulting  in 
good  engineering  estimates  for  the  residual  stresses  in  laminates. 

However,  for  an  accurate  assessment,  non-Fickian  analysis  is 
needed  and  the  formulations  presented  in  this  work  provide  a  sys¬ 
tematic  method  to  treat  this  case. 

i 
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Appendix: 

The  Analytical  Solution  to  the  Diffusion  Problems  with 
Boundary  Condition  C(±L,t)  =  C0y/tH(t) 


Employing  the  notation  of  section  3,  the  transformed  governing 
equation  for  the  diffusion  problem  is: 


dz 2  D°  ° 

(A.l) 

The  solutions  in  the  transformed  space  are 

7=7  ^  V*  cosh  qz 

°  2  V^cosh qV 

(A.2) 

M(s)  =  CoV^DtanhgI. 

sz 

(A.3) 

For  simphcity  in  future  use,  we  denote,  in  eqn.(A.2), 

.  .  cosh  qz 

^  5  s  cosh  qL  * 

(A.4) 

Generally,  if  y(s)  can  be  expressed  as 

y«  =  (a.5) 

9{s) 

where  f(s)  and  g(s )  are  polynomials  of  s,  and  if  the  degree  of  g($ )  is 
higher  than  that  of  /(s),  then  y(s)  can  be  further  expressed  as 


N 

y{s)  =  H 


n=l 


fM 

g'(a  n) 


(A.6) 


where,  an  are  zero  points  of  g(s)  and  N  is  the  degree  of  g(s). 
In  the  present  case,  f(s)  =  coshgz,  and  g(s)  —  scoshqL. 
Since 


cosh  z  — 


(A-7) 
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and  g(s)  is  indeed  one  degree  higher  than  that  of  f(s).  Furthermore, 
an  can  be  solved  from  the  equation 


s  cosh  qL  =  0 


to  give 

and 

hence 


qnL  =  ± 


ooo  =  0, 

(2  n  +  l)7rt 


D(2n  +  1)2tt2 
4jL2 

Thus  C(s )  can  be  rewritten  as: 


On  = 


(n  =  0,1,2,.-.). 


C{s)  =  Co 


y/l 

2  [^(0)  S3/2  '  ^  g'(Cn)  \A(s  “  «n)J  ' 


m  i  ,  v/(« n) 


(A.8) 

(A.9) 

(A.10) 

(A.ll) 

(A.12) 


Although  each  value  of  an  corresponds  to  two  values  (+  or  — )  of 
qn  in  equation  (A.  10),  the  choice  of  any  specific  sign  is  immaterial, 
since  either  one  would  yield  the  same  results  for  the  present  problem. 
Selecting  the  positive  sign  we  obtain 


M 

=  cosh  qz, 

(A.13) 

m 

=  coshO  =  1, 

(A.  14) 

fM 

(2  n  +  1)ttz 

=  008  2 L  ’ 

(A.15) 

A*) 

=  cosh  qL  +  -qL  sinh  qL, 

(A.16) 

AO) 

=  cosh  0  -  1, 

(A.17) 

*'(«») 

(2n  +  l)7r(  ^n+1 

(A.18) 

In  inverting  C(s)  as  expressed  in  eqn.(A.12),  we  note  that  the 
inverse  of 

7M  =  -V,!  (A. 19) 


is 


\/s(s  —  a) 

f(t)  =  -7=  e"ot  erf(Vot). 
v  a 


(A.20) 
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On  the  other  hand,  the  so  called  Dawson’s  integral  has  an  alternative 
expression  (Olver,  1974) 


F(£)  =  ~~  e  erf(£i)  =  e“*2  /  e“2<ftt. 
2i  Jo 


(A.21) 


Comparing  (A.20)  with  (A.21),  then  substituting  the  above  result  into 
eqn.(A.12)we  obtain: 


8  ^  (-l)n+1  „/(2n  +  l)7r  S  (2n  +  l)7rzl 

c^t) = coVt [1 + S (krTTFF (  2-H 008 Sr1-] • 

(A.22) 

The  total  moisture  weight  gain  can  be  readily  computed  from  the 
integration  of  the  above  expression  to  yield 

MM  =  2LCoVt  [l  -  F  '  <A-23) 
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where  departures  from  straight  lines  exceed  1%.  M(<»)  =  2L  for  the  case  of  unit  boundary 
condition. 


(%)  (i)  pm  ‘ureS  jq8i3M  amjsrow 


c 


Figure  2.  The  total  moisture  weight  gain  data  and  the  theoretical  predictions  by  no 
fitting  and  Fickian  fittings  with  two  different  saturation  levels. 


series  of  decaying  exponential  functions. 


non-Fickian  fit  and  the  "upper"  and  "lower"  Fickian  fits.  (500  hours) 


non-Fickian  fit  and  the  "upper"  and  "lower"  Ficldan  fits.  (3600  hours) 


Figure  /.  Moisture  aistnoution  across  tne  laminate  s  tnickne 
non-Fickian  fit  and  the  "upper"  and  "lower"  Fickian  fits.  (8000  hours) 
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Figure  8.  The  residual  hygro-thermal  stress  Ox  across  the  thickness  of  a  [02/902ls  laminate 
at  t  =  500  hours. 


Figure  9.  The  residual  hygro-thermal  stress  ay  across  the  thickness  o 
at  t  =  500  hours. 
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Figure  10.  The  residual  hygro-thermal  stress  Ox  across  the  thickness  of  a  [02/902]s  laminate 
at  t  =  3600  hours. 
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Figure  11.  The  residual  hygro-thermal  stress  Oy  across  the  thickness  of  a  [02/902]s  laminate 
at  t  =  3600  hours. 
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Figure  12.  The  residual  hygro-thermal  stress  Gx  across  the  thickness  of  a  [02/902]s  laminate 
at  t  =  8000  hours. 


at  t  =  8000  hours. 


